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ABSTRACT 

Context. New eclipse minimum timings of the M4.5/M4.5 binary CM Dra were obtained between the years 2000 and 
^ 2007. In combination with published timings going back to 1977, a clear non-linearity in observed-minus-calculated 

times has become apparent. Several models are applied to explain the observed timing behavior. 
Q , Aims. Revealing the processes that cause the observed O-C behavior, and testing the evidence for a third body around 

5_j . the CM Dra system. 

Methods. The O-C times of the system were fitted against several functions, representing different physical origins of 
' the timing variations. 

Results. An analysis using model-selection statistics gives about equal weight to a parabolic and to a sinusoidal fitting 
function. Attraction from a third body, either at large distance in a quasi-constant constellation across the years of 
CO ' observations or from a body on a shorter orbit generating periodicities in O-C times is the most likely source of the 

^ ' observed O-C times. The white dwarf GJ 630. IB, a proper motion companion of CM Dra, can however be rejected as 

, the responsible third body. Also, no further evidence of the short-periodic planet candidate described by Deeg et al. 

00 . (2000) is found, whereas other mechanisms, such as period changes from stellar winds or Applegate's mechanism can 

be rejected. 

04 ' Conclusions. A third body, being either a few- Jupiter-mass object with a period of 18.5±4.5 years or an object in the 

mass range of 1.5Mjup to O.IMq with periods of hundreds to thousands of years is the most likely origin of the observed 
. minimum timing behavior. 

OO , Key words. Stars: individual: CM Dra - binaries: eclipsing - Eclipses - planetary systems 

O 

> 
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CM Dra (LP 101.15, G225-067, GJ 630.1) is a detached class of planets. Motivation for this work also arises from 
■ " " ' spectroscopic eclipsing M4.5/M4.5 binary with one of the previous successes of precise timing measurements to detect 
lowest known total masses, of O.44M0. With its nearly the presence of planets. The first known extrasolar planets 
edge-on inclination of 89.59° (see iKozhevnikova et aLll2004l were detected through li ght-time effects in the s ignals of 
for the most recent orbital and physical elements) it was the Pulsar PSR1257-I-12 (jWolszczan fc Frai]||1992D and re- 
chosen as the target of the first photometric search for plan- cently, sinusoidal residuals in the pulsation frequency of the 
etary transits. Performed by the 'TEP' project, with an sdB pulsating star V391 Peg have been explained through 
intense observing campaign during t he years 1994 - 1999 the presence of a giant planet l|Silvotti et al.ll2007l) . leading 
(jPeeg et all 119981: iDovle et al] |2000[ ). over 1000 hours of to the first detection of a planet orbiting a post-red- Giant 
coverage of that system were obtained with several Im- star. In both cases, timing measurements have led to detec- 
class telescopes. This lightcurve was initially searched for tions of planets that would have been difficult or impossi- 
the presence of transits from planets in circumbinary 'P- ble to find with other planet-detection methods, a situation 
type' orbits with 5 - 60 day periods, with a negative result that is similar to the detection of planets around eclipsing 
(|Dovle et al.ll2000D . The same lightcurve provided, however, binaries - unless they exhibit transits, 
a further possibility to detect the presence of third bodies, A first a nalysis of 41 eclipse minima times for CM Dra, 
from their possible light-time effects on the binary's eclipse presented in iDeeg et all ()2000[ ) (hereafter DDKOO) gave a 

minimum times. 

^ A possible circumbinary planet is HD202206c, orbiting 
around HD202206a and b. This depends however on the clas- 
sification of the 1 7.4Miup object HD2 02206b as being a planet 



1. Introduction To date, no unambiguously circumbinary planets have 

been detecte43, and their discovery would constitute a new 
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low-confidence indication for the presence of a planet of 
1.5 - 3 Jupiter masses, with a period of 750 - 1050 days. 
This result was based on a power-spectral analysis of the 
minimum timings' residuals against a linear ephemeris, in- 
dicating a periodic signal with an amplitude of about 3 
seconds. Motivated by the result of DDKOO, we continued 
surveying the CM Dra system with occasional eclipse obser- 
vations in the following years. During this time, it became 
increasingly clear that a simple linear ephemeris would not 
provide a sufficient description of the general trend of the 
eclipse times any longer. This led to the objective of this pa- 
per - a thorough and systematic discussion of the possible 
processes acting on this interesting system, and an evalua- 
tion of the presence of a third body of planetary mass or 
heavier. 

In the following sections, we present the data for this 
analysis (Sect. [2|), evaluate the effect of several physical 
processes on the eclipse timings (Sect.|3|), and compare the 
significance of several numerical fits for the explanation of 
the observed trend (Sect. [4]). This is followed in Sect. [5] by 
a discussion of the physical implications of these fits, with 
conclusions in Sect. [6] 

2. The observational data 

The minimum times analyzed here include all minimum 
timings derived in DDKOO. Since its pubhcation, we per- 
formed dedicated eclipse observations with the IAC80 
(0.8m) and the INT 2.5m telescopes within the Canary 
Islands' Observatories and with the Kourovka 0.7m tele- 
scope of the Ural State University in Ekaterinburg, Russia. 
In summer 2004, one of the fields surveyed for plan etary 
transits by the TrES collaboration ijAlonso et al.ll2007D con- 
tained CM Dra and time series spanning s everal weeks were 
obtai ned from the Sleuth 10 cm telescope (jO'Donovan et al.l 
120041) . In all cases, photometric time series were obtained 
from whi ch the eclipse times were m easured using the 
method of iKwee fc van WoerdenI (|1956D . We accepted only 
results where the formal measurements error from that al- 
gorithm was less than 10 seconds, which required data with 
a largely complete and uninterrupted coverage of individ- 
ual eclipses. The timings obtained were then corrected to 
solar-system barycentric (BJD) times with the BARYCEN 
routine in the 'aitlib' IDL library of the University of 
Tiibingen. 

Only a few timings of CM Dra o f comparable quality 
could be found in the literature: iLacvl ([l977) gives a total of 
9 minimum times from observations in 1976. These timings 
scatter in 0-C by about 30 seconds, and were based on 
photoelectric data with frequent gaps, including several 
incomplete eclipses. We therefore chose to remeasure them 
from .Lacvi 's tabulated lightcurve with the same procedures 
used for our own data, from which only two timings of 
sufficient quality could be obtained for the further analysis. 
Only two additional timings from t he years 2004 and 2005 
could be fou nd in the literature (jSmith fc CatonI l2007t 
iDvorakI 12005) . both in good agreement with the other 
data. The final data set contains 63 minima, of which 27 
are primary and 36 are secondary eclipses. 

Since these minimum timings cover about 30 years 
of observations and are of consistently high precision, 

^ http:/ /astro. uni-tuebingen.de/software/idl/aitlib/ 
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Fig. 1. 0-C values of CM Dra eclipse minimum times from 
Table [1] against ephemcrides from DDKOO 

the effects from the 18 leap-seconds that have been 
introduced into Universal Time during that span need to 
be corrected for. Consequently, all minimum times used 
in this work were converted from the conventional UT to 
TAI (International Atomic Timej3 which is a timescale 
with constant and uniform fiow, without discontinuities 
from leap-seconds. All these minimum times arc listed in 
Table [T] together with the formal errors of the minimum 
times, the cycle number E and 0-C (ob served - cal c ulated ) 
residuals against the ephemerides of iDeeg et aD (|l998l ). 
For data newly presented in this work, the originating 
telescope is also indicated. 



3. Physical processes' effects on 0-C times 

The further analysis of the eclipse minimum times is based 
on the analysis of the temporal development of the 'O-C 
residuals between observed and calculated (expected from 
ephemeris) eclipse minimum times. In these, 

(O - C)e ^Te~ T,,e , (1) 

where T^^e refers to a minimum time calculated from an 
ephemerides at 'Epoch' or cycle number i?, and Te refers to 
the corresponding observed minimum time. The observed 
times Te are related to the minimum times in the binary 
system's rest frame by Te = T^ + cIe/c, where (Ie is the 
distance from the observer to the binary at cycle E and c 
is the velocity of light. Hence, based on a linear ephemeris 
Tc,E = EPc -\- Tcfl, where Pc and Tc,o are the ephemerides' 
period and time of conjunction. The difference, (O — C)e^ 
is given by: 

(O - C)e ^ T'e + cIe/c - BP, - r,,o (2) 

This general expression allows us now to develop the cases 
to be considered in our analysis. 

3.1. Accelerated binary systems with constant period 

First, we review the case of a binary system moving at a 
constant velocity vq relative to the observer, with vq << c. 
The distance to this system is given by ds ^ + v^EP' , 

^ The conversion TAI - UTC can be obtained from 
|ftp://maia.usno.navy.mil/ser7/tai-utc.dat 
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Sourc es: 1: Remeasured from lightcurve of lLacvl l| 19771) . 2: lDeeg et aII(|2000D . 3: This work, 4: lDvorakl (|2005D . 5: ISmith fc CatonI 
l|2007l ). 
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and the minimum times in the binary system's rest-frame 
are given by = Tq + EP' , where P' is constant. The sub- 
indices '0' refer to the parameters' values at the moment 
when E = 0. From Eq. [21 we obtain then a relation linear 
in E: 



(To 



EP' + (do + voEP')/c 
- r,,o) + EiP - P,) 



EPc - T,.o (3) 
(4) 



where P — P'{1 + vq/c) is the observable period and 
Tq =Tq + do I c is the minimum time that should have been 
observed at = 0. Hence, in a system where the observed 
0-C times deviate linearly from the ephemeris, the terms 
= (To — Tc.o) and k,\ — (P — Pc) indicate errors in the 
derivation of the original ephemerides {Tcfi,Pc), but do 
not have any further physical meaning. 

If we consider as D{E) any additional distance to the 
binary that cannot be expressed by a linear term, so that 
ds — do + voEP' + D{E), then the 0-C times will be mod- 
ified by the 'light-time effect' D{E)/c: 



{O - C)e = {To - T,^o) + E{P - P,) + D{E)/c 



(5) 



Hence, the non-linear components of (O — C)e de- 
scribe the non-linear components of the distance to 
the binary. The non-linear distance component D has 
necessarily to be caused by an acceleration process, with 

D{E) = Up a\\{t) dtdt, where a II is the acceleration 

along the line of sight to the binary. Within the scope of 
this work, the following acceleration scenarios have been 
considered: 

i) For the case of a constant acceleration a||, with 
D{E) = \a\\{EP'f we obtain now: 



(O - C)e - (To - Tc.o) + E{P - Pc) + S'aii — (6) 



= KO + Eki + E K2 



(7) 



where the Hi refer to the polynomial coefficients of the fit 
given in Tabled 

ii) In the case of an acceleration that undergoes a con- 
stant change a||, the 0-C times behave like a third order 
polynomial: 



iO-C)E = K0 + Eki + E^K2 + E^K3 

where K3 is given by: 



K3 = au 



pr3 

6c 



(8) 



(9) 



whereas the coefficients — K2 remain identical to the pre- 
vious case. 

in) For a binary that is accelerated due to third body 
on a circular orbit, the amplitude of the timing variation 
D/c is then given by: 



D/c = 



TO3 d 



||63 



ma 



{mi, + TO3) c {nib + m3)2/'^c \ 47r2 



p2Q\ 1/3 



sin i (10) 



where d||b3 is the line-of-sight-component of the distance 
between the barycenter of the binary and the third body, 
and m3 and mb are the masses of the third body and the 



binary stars, respectively. In the right hand term, d^^^s has 
been substituted for P3, the period of the third body, with 
i being the inclination of its orbital plane and G is the 
gravitation al cons t ant. A development of the general case 
is given by llrwini (Il959l). with furthe r examples of recent 
applications in iDemircan fc Buddind (|l003). 

From Eq. [TOl the 0-C times are then given by 



{0 — C)e — KO + Eki + KdCOs{EK^ + (j)o) 



(11) 



with Kd — D /c, and the term Ekij, + (po describing the phase 
of the third body, where 



H4, 



P'2tt 



(12) 



and (po is the phase of the third body at To 



3.2. Variation of the intrinsic binary period 

Here we consider the consequences of a period variation of 
a binary moving at a constant velocity. The intrinsic period 
variation dP jdE shall be small enough to be considered 
constant across the observed time span. The period at cycle 

E is then given by = J ^ dE + P^ ^ E^ + P^, and 
the times of minima in the binary's reference frame are: 



Te = / PsdE 



2^ -dE 



Tn 



EPn + Tn 



(13) 



(14) 



After converting from the times Tg to the observable times 
Te by inserting this equation into Eq. [2] , we obtain a 
quadratic equation that is similar to the accelerated sys- 
tem described by Eq. [3 with the only difference being that 
the parameter K2 is now given by: 



K2 



1^(1 + ^) 
2 dE^ c' 



(15) 



4. The observed minimum times 



As in DDKOO, 0-C times were derive d using the linear 
ephemerides given by iDeeg et al] (|1998() , which was based 
on a fit to eclipse timings observed from 1994 to 1996. For 
the present work, this ephemerides was converted to TAI 
by adding 29 seconds, which corresponds to the difference 
TAI-UT that was in effect during most of that time (July 
1, 1994 to Dec 31, 1995). The 0-C values of DDKOO differ 
therefore by a maximum of only 1 sec between the original 
work in DDKOO and the p r esent work. The ephemerides 
conversion from iDeeg et "all (|1998f) to TAI is: 

Tci = 2 449 830.75734 ± 0.000 01 -t- £^ {TAI) (16) 
Tell = 2 449 831.39037 ± 0.000 01 + Pc E {TAI) (17) 

where Td and Tdi refer to primary and secondary minima, 
respectively, E is the epoch or cycle number, and the pe- 
riod is given by Pc = 1.268 389 861±0.000 000 005 days. The 
0-C values against that ephemerides are given in Table [1] 
and shown in Fig [1] Since the writing of DDKOO, a clear 
trend of increasing 0-C values has become apparent, which 
is also apparent from the exte nsion i nto th e past through 
the inclusion of the values from lLacvl (|l977l ). which weren't 
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included in DDKOO's analysis. The linear ephemerides of 
DDKOO, therefore, no longer provides the best description 
of the observed 0-C values. However, we chose to maintain 
this ephemerides, since small errors in the parameters of a 
linear ephemerides have no physical meaning in the inter- 
pretation of higher-order 0-C dependencies, as was shown 
in Eq. [H 

4.1. Temporal evolution of the phase of the secondary eclipse 

A primary concern was assuring that the primary and sec- 
ondary eclipses were not undergoing any evolution of their 
relative phase due to, for example, variations in eccentric- 
ity or in the argument of periastron due to apsidal motion. 
Therefore we investigated if the orbital phase of the sec- 
ondary eclipse underwent any variations. This was done 
through a comparison of two sub-samples of the timings 
from Table [1] taking an early one from Epochs to 400 
and a late one from Epochs 2500 - 2900. In both samples, 
the average of the primary eclipses was set to zero, and the 
corresponding phase of the secondary eclipses were calcu- 
lated, giving: 

phase = 0.499064 ± 0.000017 for E = - 400 (18) 

phase = 0.499039 ± 0.000029 for E = 2500 - 2900 . (19) 

The difference between these two values is well within their 
error-bars; hence no relevant change in the phase of the sec- 
ondary eclipse was detected. Since primary and secondary 
eclipses of CM Dra have very similar depths, resulting in 
timing measurements of similar precision, both types of 
eclipses were treated together and equally in the further 
analysis. 

4.2. Model fits to the 0-C times 

Our initial intent was a direct fitting of the functions 
described in Section [3] to the 0-C residuals of Table [1] 
However, in the subsequent statistical analysis we noted 
two problems: First, the average formal error in the indi- 
vidual 0-C measurements from Table [T] is 3.90 seconds. On 
the other hand, the standard deviation among the residuals 
(observations - fits) could not be reduced below 5.6 seconds. 
Reduced values were not less than Xrod ~ 4.2, even when 
applying 5th or 6th order polynomial fits, whereas a 'good' 
fit should indicate values of Xrcd ~ 1- Since the fits are 
not intended to - and cannot - model the point-to-point 
variations and furthermore, since the presence of periodic 
short-frequent 0-C timing variations can be ruled out (see 
section [S75|) ■ we concluded that the formal errors are sub- 
estimating the real measurement errors. In the further anal- 
ysis, a value of 5.6 seconds was adopted as a minimum mea- 
surement error and errors smaller than this were set to 5.6 
sec. With fourth to sixth order polynomial fits to these data 
we now obtain Xtcd ~ 1- 

A second problem arose because the data-points aren't 
uniformly or randomly distributed, but clustered into 
yearly observing seasons. These clusters act like pivots for 
the fitting functions, reducing the degrees of freedom of the 
model fits over the number of data points. The choice of 
the correct number of degrees of freedom is, however, im- 
portant for the model comparison performed in Sect. 14.41 
Consequently, for each year of observations, we generated 
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Fig. 2. Seasonally averaged 0-C values of CM Dra, with 
parabola fit (dashed line) and sine-linear fit (solid line) 



one single data point from a weighted average of each sea- 
son's points. The resulting binned 0-C times are shown in 
Fig. [2 

Fits of the functions that have been discussed in Sect. [3] 
were then performed against this sample, with their best-fit 
parameters given in Table O The fits were obtained with 
the IDL library routines "POLY JIT" for the pol ynomial 
fits, and "AMOEBA" (based onl Press et all |l992^ for the 
sinusoidal fit, with both algorithms performing minimiza- 
tions of the x^ values. It should be noted that the aim of 
the sine-linear fit was to test if the data can be well fit to 
one or a few cycles of such a function, with a period longer 
than about 10 years. For a search for shorter periodicities 
see Sect. 15.31 

The parameters' errors given in Tableware 1-sigma con- 
fidence limits, whose derivation is described in the next sec- 
tion. Furthermore, the Table's last three columns give the 
best- fit values, calculated for a sample standard devia- 
tion of s = 3.2 s, the Akaike Information Coefficient AICc, 
and the Akaike weights Wi ; both of which will be introduced 
in Sect. SH 

4.3. Errors of fit-parameters 

For an estimation of the errors of the fit-parameters we ap- 
plied initiall;Y_t|i6_comn^ resampling or bootstrap method 
fe.g. lCameron et al]|2007f ). consisting of repeated model fits 
to a synthetic data set. It requires one to assume some ref- 
erence function that describes correctly the general trend 
of the data, against which residuals of the data points are 
generated. The synthetic data are generated by permutat- 
ing the residuals among the data points. The model to be 
evaluated is then fitted against the synthetic data and dis- 
tributions of the obtained fit-parameters are used to es- 
timate the likelihood distribution of the parameters from 
the model-fit on the original data. This method is easy to 
implement and leads to synthetic samples with properties 
similar to the original data, but in this work's context two 
problems arose: 

First, the analysis is based on the assumption that the dis- 
tribution of the fit-parameters on the synthetic sets is iden- 
tical to the likelihood distribution of the parameters ob- 
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Table 2. Parameters of fits to 0-C times 





s 


lO^^s 


lO^^S 


^3 


S 


^0 <pQ 

lO^^rad rad 


AIC 




Linear 

O - C = K0 + Eki 


3.6 ± 1.0 


7.3 ± 0.4 








37.94 


49.4 


0.002 


Parabola 

O - C = Kq + Eki + E^K2 


-1.5 ± 1.5 


8.7 ± 0.5 


7.0 ± 1.3 






11.57 


37.4 


0.649 


Third Order Polynomial 

O - C = Kg + Eki + E^K2 + E^K2 


-1.1 ± 1.7 


8.0 ± 2.0 


7.8 ± 2.6 


0.4 ± 1.1 




11.45 


41.6 


0.079 


Linear + Sine 

O — C — kq + Eni + s\n{EK^ + 4>q) 


7 3+2.1 

' •■^-4.2 


5 i 






°- ' -3.0 


l.lSt°il i-^-Vr 6-18 


39.1 


0.270 



tained from the fit of the original data, which is far from 
certain. 

Second, resampling is based on the assumption that the ref- 
erence function is a correct description of the trend of the 
data, and that residuals against it are measurement errors. 
This approach may be justified if the underlying physical 
model - and the function that describes it - is known, and 
only a refining of parameters is required. In this work how- 
ever, we are also faced with an uncertainty about the nature 
of the reference function. 

A method that overcomes these problems, and which 
has come to the awareness of astrophysicists in recent 
year s, is the Markov Chai n Monte-Carlo (MCMC) method 
Ce.g. iTegmark et al.|[200l: iFordI [20051 : iHolman et all 120061: 
iBurke et al.l l2007f ). MCMC is based on the states of a 
Markov Chain undergoing random variations, whose prob- 
abilities are however directed by the maximum likelihood 
estimator (MLE) of the fit-parameters (relative to the orig- 
inal data) at each step in the chain. The frequency of states 
of the Markov Chain at a given parameter value indicates 
then the posterio r pro b abilit y distribution of that param- 
eter. We refer to iFordI (|2005[ ) for further references to the 
MCMC, as well as for the implementation of the MCMC 
with the Metropolis-Hastings Algorithm that was used in 
our data analysis. 

Our final error- analysis is however not based on the 
posterior probability distributions derived from histograms 
of the parameter distributions. Instead we have used the 
MCMC as a tool to explore the re lation of against the 
multivariate parameters. Similar to lBurke et al.l ( 20071 ). for 
any recorded step of the MCMC chain, the encountered pa- 
rameters and the corresponding values were registered. 
As an example. Fig. [3] shows the relation between the 
parameter in the sine-linear model fit against x^- There is 
a clear lower limit of for a given parameter value. With 
increasing numbers of iterations in the Markov Chain, this 
lower limit approaches the best possible fit at a given pa- 
rameter value. Hence, the lower limits of may be used as 
tracings of the best-fit against one or multiple param- 
eters. Since the distribution-density of a MCMC sequence 
gravitates towards the regions of lowest x^, the MCMC 
method may therefore be used as a simple tool to trace 
low-sigma confidence regions around the best-fit parame- 
ters. This application of the MCMC also avoids a problem 
that easily occurs in the interpretation of posterior prob- 
abilities from the parameter densities. As can be seen in 
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0.0000 



0.0005 



0,0010 



0.0015 



0,0020 



Fig. 3. Distribution of the parameter from the sine- 
linear model against x^ from a Markov Chain of 20 million 
steps. The dashed lines indicate the parameter's 1-ct confi- 
dence region, corresponding to Ax^ — 1 over the best fit's 
X^ value (cross) 



Fig. [21 the distribution of points also has dense zones close 
to the left cutoff, although the minimum value of x^ in 
that zone corresponds to very poor fits. With the flat de- 
pendency of the best-fit x^ against in that zone, the 
Markov Chains had difficulties returning to better-fitting 
values. In the given case, the Markov Chains instead went 
into 'exploring' a wide multi-parameter space of poor mod- 
els that opened up close to the left cutoff in Fig. [31 

For the errors indicated in Table [21 we used chains with 
a length of 20 million iterations, after discarding the first 
1000 steps of the chains. 



4.4. Model selection 

After applying the fits of several models, physical interpre- 
tations will need some information about the likelihood that 
a given fit corresponds to the true behavior of the observed 
data. This question is commonly referred to as 'model se- 
lection'. Values such as x^ (see Table [2]), or the 'reduced 
Chi-square' of x^ jv., where v is the degrees of freedoms. 
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give a general indication of the quality of a fit. However, 
in the case of small or moderate differences of fit-quality 
among models, they serve little to generate statements that 
are useful for the model selection. 

Probably the best-known method of assigning likeli- 
hoods in fit-comparisons is the 'F-test'. It is however valid 
only for the comparison of nested models, such as poly- 
nomials of different orders, and was therefore n ot used 
further. The Akaike Information Coefficient (AIC, lAkaikd 
Il974[ ) does not have this limitation, which led to its use in 
the further investigation. For an introduction to the use of 
th e AIC as a tool for model sel e ction, we refer the reader 
to iBurnham fc Anderso^ (|2004f ): iLiddld (|2007[ ): iMazerolld 
(|2004l ). Using the residuals' squared sum (RSS) as the like- 
lihood estimator, where RSS = J^iVi ~ fi)^ with yi being 
the data and /j the model values, the AIC is given by: 



AIC = -nln(RSS/n) + 2A: 



(20) 



where n is the number of observations and k is the num- 
ber of model parameters. Here we use the generally pre- 
ferred (jBurnham fc Andersonl[2b04[ ) second order corrected 
coefficient 'AlCg', which is valid for both small and large 
samples: 



AIC^ = AIC 



2k{k + l) 



N 



1 



(21) 



While smaller AIC values indicate better fits, their abso- 
lute values don't have any meaning. They are only useful 
if values from different models are compared. If the best 
among several models has a value of AIC cmim then for 
any model i the differences = AICc,i — AlCc.min may 
be calculated. Differences of < 2 indicate substantial 
support (evidence) for model i; models where 4 < A^ < 7 
have considerably less suppo rt, while models with A^ > 1 
have essentially no support (jBurnham fc AndersonI [200^ . 
For the comparison within a set of R models, normalized 
Akaike weights may be derived for each model i, with 



exp(-A,/2) 
Etiexp(-A./2) 



(22) 



where all weights Wi sum up to 1 (see Table [2]). Following 
lAkaikd (1981), these weights may be interpreted as a likeli- 
hood that can be assigned to each of the models, with the 
parabolic fit being most likely, followed closely by the sine- 
linear fit. A word of cauti on, also refiec ted in several refer- 
ences about this topic (e.g. Liddle'2007), should however be 
given against the use of these statistical values as a strong 
argument in favor of one or the other model: There are sev- 
eral alternative indicators available, such as the Bayesian 
Information criterium (BIC, Schwarz 1978 ) or the De viance 
Information Criterion (DIC. lSpiegelhalter et al.ll2002[ ). with 
different 'penalizations' for models with additional degrees 
of freedom. While these criteria indicate similar preferences 
for models with well separated or AICc values, inter- 
changes in ranking may happen among models that are 
close in AICc values. This cautionary position was backed 
by a calculation of the BIC for our models. In that case, the 
sine- linear fit was ranked best, with a slightly lower (and 
hence 'better') BIC than the parabola fit, while the ranking 
of the other models remaining unchanged. 

In summary, both the simple linear fit and the third 
order polynomial fit have significantly less support than 



the parabolic fit and the sine- linear fit. In Section [5] we 
therefore focus on the physical implications from these two 
top-ranked models. 

5. Discussion: Possible causes of the observed 0-C 
times 

Common to all fits, the linear parameter ki has fairly simi- 
lar values indicating clearly that the average period of CM 
Dra across the recorded o bservations i s seve ral millisec- 
onds longer than given by iDeeg et all (|1998D . As shown 
in Sections 13. II and [3. 21 a parabolic 0-C function may arise 
from two causes, an intrinsic variation in the system's pe- 
riod, or a light-time effect from a constant acceleration of 
the entire binary system. In both cases, the only interesting 
parameter is the quadratic term, found to be (see Tabled 
K2 = (7.0 ± 1.3)10~'^ sec/period. 

5.1. Intrinsic period variation 

Considering an intrinsic period variation, the change in 
period- length per cycle is given by EqlTS] with ^ << 1 
as: 

dP' 
'dE 



2k2 = (1.4 ± 0.26) X lO^'^s/period, 



(23) 



The corresponding unitless period change per time is given 
by: 



dP _dP 1 

~9r ~ aE 



1.28 X 10"" 



(24) 



iDemircan et al.l (|2006D performed a statistical study of 
the orbital parameters of a sample of detached chromo- 
spherically active eclipsing binaries of different ages. Out 
of that sample, they concluded that their periods decrease 
with an average value of a = 3.96 x 10~^°yr"^, with a de- 
fined by the differential equation dP/dt = —aP. This value 
may be considered constant throughout a large part of a bi- 
nary's evolution and is due to angular momentum loss from 
magnetically driven stellar winds. The corresponding pe- 
riod change of CM Dra would be ^ = -1.38 x lO'^^^ ob- 
tained by multiplying —a with CM Dra's period in units of 
years (3.47x 10~'^yr). This value is of opposite sign than the 
observed one (Eg. [24)1 . and is an order of magnitude smaller. 
Hence, angular momentum loss may well be present, but is 
not detectable in the current minima timings. 

5.2. Acceleration due to a quasi-stationary attractor 

The second source for a parabolic shape in an 0-C dia- 
gram could be a constant acceleration of the binary along 
the line of sight, with changes of acceleration strength and 
direction during the span of observations being negligible. 
The acceleration term is given from Eq. [7]as: 



2^ ^ (3.5 ± 0.7) x 10" V/s^ 
P ^ 



= (1.17±0.22) X 10"''s 



17„-1 



(25) 



(26) 



We note that this acceleration is at least two or- 
ders of magnitude larger than the acceleration of the 
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Solar system, currently constr ained within a fewx lO^^^s^^ 
(jZakamska fc Tremaind [20051) . A quasi-constant accelera- 
tion may be caused by a third body of mass at a dis- 
tance far away enough so that mutual orbital motions don't 
lead to significant changes in the acceleration vector. The 
acceleration on the binary caused by a third body is given 
by: 



Cms 



(27) 



where r is a distance vector from the barycenter of the bi- 
nary towards the third body. We note that this acceleration 
is independent of the mass of the binary. The acceleration 
component along the line of sight is then given by 



Cms 



■ cos I sm I 



(28) 



where i is the inclination, defined here as the angle between 
r and the plane of the sky, and r± = r cos i is the lateral 
component of r. The equation above gives ay =0 for in- 
clinations of both 0° and 90° and a maximum for a|| at 

inclinations of i — arctan = 35.26°, leading to 



< 0.3849 



(29) 



Converting to common astronomical units, the minimum 
mass 7713 at a given lateral distance is then obtained by 



m3 
M0 



> 438.26 



(^) 



2^ (au) 



(30) 



Replacing r± by an angula r separation based on th e 
distance to CM Dra (15.93 pc; IChabrier fc Baraffdll995[ l. 
and using a|| = 3.5±0.7x 10~*m/s^, we may now derive the 
minimum mass of a possible third body at a given angular 
separation from CM Dra: 



> 0.0030 f— 

V arcsec/ 



(31) 



The third order polynomial fit, while resulting in a lower 
weight in the model selection, is not to be discarded com- 
pletely. It would describe a system undergoing a constant 
change in acceleration a|| . There is however no physical pro- 
cess that generates a truly constantly varying acceleration. 
Hence the third-order polynomial fit may only describe 
cases where a|| is quasi-constant across the observing time- 
span. The third order polynomial fit may, however, provide 
the first terms of a Taylor-expansion of the true acceleration 
process of the system. For a slowly changing acceleration, 
like a cyclic one with long periods of 0(100 yr) or longer, 
the 3rd order polynomial may therefore give a better de- 
scription of the 0-C times than the parabolic fit does. From 
our fit, however, with the value of K2 being very close to 
the one from the parabolic fit, the derived acceleration a|| 
and the constraint for a third mass given in Eq. [2l] do not 
significantly differ. 

A possible source for an acceleration of the CM Dra 
system may be the nearby white dwarf GJ 630. IB 
(WD 1633-f57, LP 101.16, G225-068). This has long been 
recognized a. s a p roper-motion companion to CM Dra 
(|Giclas et al.l llQTll ) at an angular distance of 26", which 



corresponds to a lateral distance r± of 414AU. With a pe- 
riod of O(IO^) yr, the criteria of a quasi-constant acceler- 
ation vector during the 31 years of observational coverage 
is clearly given. Following Eq. I31[ a minimum mass of ms 
of 2.0 ± O.4M0 is however obtained for the white dwarf 
- much above the typical white dwarf masses of 0.5 - 0.7 
Mq, and clearly above the Chandrasekar limit for white 
dwarfs of « 1.4Mq. Assuming a mass of 0.6AfQ for GJ 
630. IB, this object would contribute an acceleration of only 
a ^ 1 X 10~^m/s^ on the CM Dra system. 

While GJ 630. IB can be ruled out as a source of the 
observed 0-C variations, they may be caused by still undis- 
covered bodies in the brown-dwarf mass regime (13 to 80 
-^^jup) at maximum distances of about 5 arcsec, or by a 
planetary-mass object (with less then 13 Jupiter masses) 
at a maximum distance of 2 arcsec. 



5.3. An orbiting third body 



As shown in section 14.41 the sinusoidal fit matches the 
observed 0-C times about as well as the parabolic one. 
From the fitted parameter, k^, we obtain with Eq . [T^ and 
ma << mt, a period of 



Pa = = 18.5 ± 4.5yr 

K0 



(32) 



Eg. llOi with D/c— Kd can be rewritten as: 

2.1Mj,p . (33) 



ma sm i — Kd 



f nib / Ps^ 



\Mq I yr 



(34) 



With the above value for Pa and the fitted one for k^, this 
leads to: 



ma sin I = 1.5 ± 0.5Afjup 



(35) 



We note that such an object would have an orbital half-axis 
of about 5.3AU, with a maximum separation from CM Dra 
of about 0.35 arcsec. 

While the sinusoidal fit indicates a periodicity on a time- 
scale of 20 years, we also performed a search for higher 
frequencies. For this, the same sine-wave fitting algorithm 
used in DDKOO was employed. This analysis was performed 
on the 0-C residuals using the parabola fit, and included 
only the relatively dense surveying that started in 1994, 
thereby excluding the two isolated early values obtained 
from Lacy (1977). The resulting power-spectrum is shown 
in Fig. m 

The single broad peak at a period of 950 da ys and with 
an amplitude of 3 seconds that was apparent in iDeeg et "al] 
(|1998 Vs Fig 2.C has now disappeared, and been replaced by 
several peaks with amplitudes close to 3.5 seconds. We note 
that this amplitude is close to the sample standard devia- 
tion of the yearly averaged data of s = 3.2s. Since no single 
peak is outstanding, no indications for any periodicites on 
timescales of ^10 years remain. 

5.4. Applegate's mechanism 

We also evaluated t he m echanis m introdu c ed b y 
lApplegate fc PattersonI (|1987[ ) and lApplegatd (|1992[ ). 
which may give rise to orbital period modulations of 
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Fig. 4. Power-spectrum obtained from residuals of 0-C val- 
ues against the parabola fit 



binaries from the periodic variation of the shape - and of 
the quadrupole moment - of a magnetically active binary 
component across its activity cycles. The notion that a 
component of CM Dra may be magnetically active cannot 
be completely discarded, since sev eral large f l are events 
have been reported in the literature (jLacvlllQTTi : iKim et aTl 
fl997t iDeeg etallllQQl . 

In Applegate's model, the angular momentum of the en- 
tire system remains constant, but the distribution of angu- 
lar momentum between the components' mutual orbit and 
the active star's internal rotation varies. Energy taken up 
by the variation in the internal rotation has to be reflected 
in the active star through a corresponding luminosity vari- 
ation. Applegate's original calculation considers the energy 
taken up by differential rotation between an inner stellar 
core and an outer shell of O.IMq, something which is in- 
appropriate for either com ponent of CM Dra, wit h masses 
of 0.207 and 0.237 Mq (|Kozhevnikova et al.|[200l . respec- 
tively. We foll owed therefore the more general calculation 
introduced bv iBrinkworth et "aTl (|2006D . which leads to the 
rotational energy that has to be provided in order to explain 
a given period change, while doing this for any distributions 
of the stellar mass into core and shell. The period variation 
under consideration is given by (iApplegat(j|l992f) : 



(36) 



where Kd is the amplitude of the 0-C variation, P the bi- 
nary period and P3 the modulation period. With corre- 
sponding values taken from Table [5] and Eq. [321 we obtain 
AP — 0.010 seconds. The minimum energy to produce such 
a period change for any distribution of stellar mass into 
core and shell (assuming a mass distribution following the 
Lane-Emden equation for a polytrope of n=1.5) amounts to 
3.8 X 10"'^ erg if CM Dra A is considered as the active star 
(see Fig. [5]); it would be shghtly higher (4.4 x 10''^ erg) in the 
case of CM Dra B. This energy may be compared to the lu- 
mino sity of CM Dra of 1.06 x 10~^Lq (jKozhevnikova et al.l 
l2004f ). which corresponds to a release of radiant energy of 
2.4 X 10'"' erg over the same span of P3 = 18.5 years. With 
the energy required for the period change being two orders 



Fig. 5. Energy required to vary the internal rotation of CM 
Dra A in order to reproduce the observed period change 
with Applegate's model, for all possible values of CM Dra 
A's shell mass. The lowest amount of energy required is 
3.8 X lO'*^ erg at a shell-mass of O.OI6M0 



of magnitude larger than the radiant energy, Applegate's 
mechanism can definitively be discarded as a source of the 
period variations. 



6. Conclusions 

The 0-C timing shows a clear indication of a non-linear 
trend. After a review of potential causes (previous section), 
the remaining explanations are given by the presence of 
a third body. The nearly equal statistical weight of the 
parabolic and the sinusoidal fits currently prevents setting 
a clear preference for the one or the other model. For the 
period of the third body there exist two distinct possibili- 
ties: A Jupiter-type planet with Msini of 1-2 Mjup with a 
period of 18.5 ± 4.5 yr, or an object such as a giant planet 
or heavier, with a period of hundreds to thousands of years. 
Intermediate-length periods of about 25 - 100 years are less 
likely due to poor compatibility with either the sinusoidal or 
the polynomial fits. Regarding shorter periodic bodies, the 
power-spectrum shown in Section [5.31 excludes 0-C modu- 
lations with amplitudes of larger than « 3.5 seconds with 
periods of ^10 years. The sensitivity of 0-C timing de- 
tections against orbiting third bodies decreases with their 
period, and hence Jupiter-mass objects on such shorter pe- 
riods cannot be excluded. A low-significance candidate for 
such an object with a period of about 900 days was pre- 
sented in DDKOO. That candidate was based on a single 
peak in a power-spectrum with an amplitude of 3 seconds, 
whereas the newer data show several peaks of amplitudes 
up to 3.5 seconds. The newer data therefore cannot refute 
that candidate, however it has become most likely to have 
been the result of a fortuitous combination of 0-C timing 
values. 

For long-period bodies at a quasi-constant distance and 
position during the observed time-span, Eq. [3T] allows us 
to set a maximum lateral separation from CM Dra for any 
given third-body mass. We also assume that any nearby 
body larger than about « O.IMq would have become ap- 
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parent in existing images of the CM Dra field, for which a 
larg e collection of CCD im ages exist from the TEP project 
(fPee g et al. 19981: lDovle~e t al. 2000), or in 2Mass images in 
the IR. A maximum lateral distance of 6 arcsec, or 95AU 
may therefore be set for the presence of an as yet undiscov- 
ered third body of ^ O.IMq. The corresponding maximum 
distances for undiscovered brown dwarfs or planets that 
could explain CM Dra's 0-C timing behavior are 5 and 2 
arcsec, respectively. While the setting of third-body mini- 
mum masses (and implied brightnesses) for a given lateral 
distance aids in the definition and interpretation of observ- 
ing projects, we note however that the observed acceler- 
ation term may be caused by relatively small third-body 
masses. If we take 100 years as the shortest circular period 
that may mimic the quasi-stationary case, such an object's 
orbital distance would be 16.4 AU. If it is aligned such 
that |r| « r|| (e.g. i close to 90deg), then solving Eg. \Tl\ 
indicates a mass of about l.SMjup. This mass may be con- 
sidered the absolute minimum mass for a third body at a 
quasi-stationary distance, with objects at larger distances 
requiring larger masses. 

In conclusion, two possibilities for the source of CM 
Dra's timing variations remain valid: A mass of a few 
Jupiters on a two decade-long orbit, or an object on a 
century-to-millenium long orbit, with masses between 1.5 
Jupiters and that of a very low mass star. Continued obser- 
vations of the timing of CM Dra's eclipses over the next 5 - 
10 years should, however, be decisive regarding the contin- 
ued viability of the sinusoidal- fit model, and hence, about 
the validity of a jovian-type planet in a circumbinary orbit 
around the CM Dra system. 
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